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ABSTRACT 

Aims. We present a stand-alone software (named EXSdetect) for the detection of extended sources in X-ray images. Our goal is to 
provide a flexible tool capable of detecting extended sources down to the lowest flux levels attainable within instrumental limitations, 
while maintaining robust photometry, high completeness, and low contamination, regardless of source morphology. EXSdetect was 
developed mainly to exploit the ever-increasing wealth of archival X-ray data, but is also ideally suited to explore the scientific 
capabilities of future X-ray facilities, with a strong focus on investigations of distant groups and clusters of galaxies. 
Methods. EXSdetect combines a fast Voronoi tessellation code with a friends-of-friends algorithm and an automated deblending 
procedure. The values of key parameters are matched to fundamental telescope properties such as angular resolution and instrumental 
background. In addition, the software is designed to permit extensive tests of its performance via simulations of a wide range of 
observational scenarios. 

Results. We applied EXSdetect to simulated data fields modeled to realistically represent the Swift X-ray Cluster Survey (SXCS ), 
which is based on archival data obtained by the X-ray telescope onboard the Swift satellite. We achieve more than 90% completeness 
for extended sources comprising at least 80 photons in the 0.5-2 keV band, a limit that corresponds to 10~ 14 erg cm -2 s _1 for the 
deepest SXCS fields. This detection limit is comparable to the one attained by the most sensitive cluster surveys conducted with much 
larger X-ray telescopes. While evaluating the performance of EXSdetect, we also explored the impact of improved angular resolution 
and discuss the ideal properties of the next generation of X-ray survey missions. 

Key words, surveys; cosmology: observations; X-rays: galaxies: clusters; galaxies: clusters: general; methods: statistical; techniques: 
image processing 
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t^j- I 1 . Introduction phenomena has been opened to direct investigation in the field 

_ i ■ , , , . of galaxy groups and clusters. The most noticeable discoveries 
Cj . In the past decade, X-ray astronomy reached an unprecedented indude the interaction between the ICM and the re i ativistic jets 
r-H . level in imaging and spectroscopic performances mostly thanks of the central radio j wMch is directl observed as bubbles 
CN to the Chandra and XMM-Newton telescopes. Their high imag- inflated by the jets; the presence of cold fronts; and the physics of 
p. ing quality and large effective area al owed us to discover an cool . cores and their metal distribut i on . However, less attention 
> increasing complexity in the morphology of X-ray extended has been id tQ the detection of faint extended sources, 
sources, which reflects the rich physics involved. Extended X- d ite the increasi interests in distant and clusters of 
ray emission may be associated with a variety of Galactic ob- laxies Based Qn the serendipitous d i scover ies of high-z mas- 
g jects: supernova remnants, star-forming regions, or planetary siye dusters of laxies th h x observation of XMM- 
nebulae. In extragalactic fields, it is mostly associated with the NewtQn and Chmdm or th h the Sunyaev-Zeldovich effect 
hot intra cluster medium (ICM) that permeates the potential well (see Jee et al 2009; Fol et al 2Qn . Santos et al 2Q1 1}> we 
of galaxy clusters and groups. Other extragalactic, extended X- t a si ificant population of massive clusters at very high 
ray sources are given by the inverse Compton scattering from redshift (z „ L5)> whose nce was considere d unlikely just a 
reiativistic jets in radio galaxies, from supernova remnants, hot fcw The detection of a sizable number of hi h . z mas . 
gas and massive X-ray binaries m star-forming galaxies, and siye dusters ye the tQ a substantial revision of the 
low-mass X-ray binaries in normal galaxies. In this paper, we standard ACDM cosmological model (Harrison & Coles 2012; 
mainly focus on galaxy clusters and groups. Hoyle et al 2012; Waizmann et aL 2012) . Moreover, a largely 
Thanks to the development of tools for analyzing the mor- unexplore d population of small or medium-z groups awaits sys- 
phology of extended sources and the capability of simultaneous tematical study . Overall, a thorough investigation of the low- 
imaging and spectroscopy with CCD, a wealth of new physical flux> extended X-ray sources is expected to provide significant 
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Robust detection and characterization of faint extended X- 
ray sources is very difficult. In the photon- starving regime of 
X-ray images (Hobartetal. 2005), Poissonian fluctuations in 
the sparse background can produce spurious source detections 
at low fluxes, which makes the firm detection of faint sources 
very difficult, irrespective of their extent. In addition, detect- 
ing an extended source is much more difficult than detecting an 
unresolved source with the same number of net photons. First, 
extended emission is spread over more pixels, resulting in a 
much lower contrast to the background, and can be more eas- 
ily swamped by the sparse background. Second, faint extended 
sources can be confused with faint unresolved sources. 

In addition, the mere detection of a source is not suffi- 
cient, since a robust characterization of its extension is also 
required on the basis of X-ray data. In principle, the only re- 
quirement to characterize a source as extended is that its mea- 
sured size must be larger than the instrument point spread func- 
tion (PSF). However, given the large number density of active 
galactic nuclei (AGNs) and distant star-forming galaxies, unre- 
solved sources largely outnumber extended ones, especially at 
low fluxes. The risk of finding faint extended emission encom- 
passing several unrelated unresolved sources is high. Moreover, 
in the classical Wolter-type design of X-ray mirrors, the PSF 
varies significantly across the field of view (FOV), introducing 
more uncertainties. Finally, extended X-ray sources are complex 
in morphology, because they come in a wide variety of shapes 
and surface brightness distributions, with scales ranging from 
a few arcsec to a few arcmin. Therefore, a considerable effort 
should be made to search and characterize extended sources in 
X-ray images. 

In this paper, we aim at providing a flexible and efficient 
detection algorithm to identify extended X-ray sources down 
to a low flux level to exploit the large amount of data in the 
archive of current X-ray missions (Chandra, XMM-Newton, 
Swift, Suzaku), and to explore the scientific cases of future X- 
ray facilities. We mostly focus on the detection of diffuse extra- 
galactic sources, namely groups and clusters of galaxies. 

The identification of extended sources consists of two 
parts: source detection and source characterization. For the 
first part, we use the same algorithm - a combination of 
Voronoi tessellation (VT) and friends-of-friends (FOF) - as 
used in the classical software vtpdetect based on the work 
by Ebeling & Wiedenmann (1993). VT is a useful tool to 
deal with typical X-ray images, which are largely dominated 
by empty pixels. A recent application has been presented by 
Diehl & Statler (2006), who used weighted VT in adaptive bin- 
ning of X-ray images. This software is not designed for source 
detection, however. 

Although it has been widely used for detect- 
ing galaxy clusters in optical images (for instance in 
Ramellaetal. (2001), Kimetal. (2002), Panko & Flin (2005), 
van Breukelen & Clewley (2009), and Barkhouse et al. (2006)) 
\ the method of VT+FOF is seldom used in X-ray cluster 
surveys. Presently, the only X-ray source detection software 
based on VT available to the community is vtpdetect 
as a part of the ciao 2 software developed for the Chandra 
mission. However, so far it has been applied mainly to ROSAT 
data (Scharf etal. 1997; Ebeling etal. 1998, 2000), and to our 
knowledge, only once to Chandra data (Boschin 2002). The 



1 These works identify overdensities in the field of galaxies using 
a third information in addition to the 2D galaxy coordinates, such as 
redshift, magnitude, or color. 

2 http ://cxc . harvard. edu/ciao4 . 4/ 



detection algorithms commonly used in X-ray surveys are 
based on wavelet transform methods (Rosati et al. 1998, the 
ROSAT Deep Cluster survey ; Vikhlinin et al. 1998, the ROSAT 
160d survey;Burenin et al. 2007, the ROSAT 400d survey; 
Fassbender et al. 2011, the XMM-Newton Distant Cluster 
Project; Lloyd-Davies et al. 2011, the XMM Cluster Survey; 
Pacaud et al. 2006, The XMM Large-Scale Structure survey; 
Barkhouse et al. 2006, the Chandra Multi-wavelength Project). 

Very few alternative lines of research are being developed for 
exploiting the entire information contained in X-ray images. So 
far, the only method fully employing the information of the pho- 
ton counts per pixel and of the counts in neighboring pixels is 
the background/source separation technique (Guglielmetti et al. 
2009). This method uses the Bayesian probability theory com- 
bined with a two-component mixture model. This way, back- 
ground and source intensity can be estimated jointly. First re- 
sults based on this technique and applied to deep Chandra fields 
are currently under scrutiny (Guglielmetti et al. in preparation). 
Given this framework, it is well worth putting a significant effort 
into developing diversified detection algorithms for extended 
sources. 

In this paper, we further develop the VT+FOF algorithm, 
making a great effort against its major shortcomings. We provide 
the flexible and efficient Extended X-ray Source detection soft- 
ware (EXSdetect), which can be easily applied to a wide range of 
X-ray images, and is optimized to detect extended sources. The 
software, written in Python and made available to the community 
on a public website, can be used for a double purpose: exploit- 
ing the rich (and increasing) archive of current X-ray missions; 
and investigating scientific cases of the next-generation X-ray fa- 
cilities, in particular those that will perform large-area surveys. 
Here we do not consider supplementary information, although 
we may use the CCD spectral information to distinguish thermal 
(X-ray soft) from nonthermal (X-ray hard) emission. This option 
is not practical, because mostly the signal-to-noise ratio (S/N) is 
very low in the hard band (above 2 keV) where the difference be- 
tween thermal and nonthermal spectra is highest. We also ignore 
the possibility of cross-correlating X-ray images with images in 
other wavebands like optical and IR. This process is very effec- 
tive in finding distant cluster candidates, but it would severely 
undermine the possibility of clearly defining the surveyed vol- 
ume, an essential requirement to derive the physical density of 
the cluster population as a function of the cosmic epoch. Our 
goal is to obtain samples selected entirely by their X-ray prop- 
erties and characterized by clearly defined completeness crite- 
ria. This requirement is mostly relevant for cosmological tests, 
where we need a robust estimate of the search volume of the 
survey to measure the comoving number density of sources. 

The paper is organized as follows. In §2 we review the 
vtpdetect algorithm and describe the improvements and the 
additions of EXSdetect. In §3 we discuss the practical case of 
the ongoing Swift X-ray Cluster Survey (SXCS, Tundo et al. 
2012), presenting extensive simulations aimed at testing the per- 
formance and evaluating the completeness of the survey. In §4 
we briefly discuss the effect of improving the angular resolution. 
Finally, in §5 we summarize our findings. 

2. Algorithm 

2. 1 . VT+FOF as implemented in vtpdetect 

There are several ways to find photon density enhancements, and 
therefore identify source candidates in an X-ray image. The sim- 
plest way is to search across the image inside sliding boxes with 
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different sizes and shapes. A more effective option is to search 
for overdensities on different scales in Fourier space, making 
use of wavelet transform. Another possibility is to calculate the 
local flux density across the image in Voronoi cells and locate 
excess in the flux density distribution. This method has been in- 
troduced by Ebeling & Wiedenmann (1993) and is implemented 
in the ciao task vtpdetect . This algorithm is a combination 
of the classical VT and FOR The classical VT is a nonparamet- 
ric method that uniquely partitions the infinite plane into con- 
vex polygon cells according to a list of starting locations (named 
Voronoi sites) in the plane. Each cell contains only one Voronoi 
site and includes all points that are closer to this site than to 
any other. The positions of all occupied pixels are assumed as 
Voronoi sites, and the local photon density value associated to 
each occupied pixel is computed as the number of photons in 
the occupied pixel divided by the Voronoi cell area. In this way, 
the raw image is converted into the local photon density map, 
in which high-density regions show up by a direct comparison 
between the photon density map and a background map. Then, 
the FOF method is used to identify source candidates. 

The VT+FOF method has two immediate advantages. First, 
it is applied to the X-ray photon event list or to the original un- 
binned image, and it preserves the full angular resolution at each 
step, at variance with other methods, where image-binning is of- 
ten adopted to avoid oversampling of the PSF and improve the 
low count statistics. Second, the VT+FOF algorithm does not re- 
quire assumptions about the shape and size of the structure one is 
looking for. In the simplest approach, a measure of the size (area) 
of the source candidate and the corresponding aperture photom- 
etry is directly obtained by defining the source edges as the loci 
where the photon density value is equal to the background den- 
sity, without using a predefined aperture region. 

Despite these properties, vtpdetect is not widely used 
compared to wavelet-based algorithms. The construction of 
VT has been considered to be very time-consuming, and 
vtpdetect was suggested to be applied only to small fields 
or low-event-density regions (Ebeling & Wiedenmann 1993). 
Another major reason for this infrequent use is that the FOF al- 
gorithm has a severe source-blending problem, i.e., the merging 
of two or more neighboring sources into one. This is a well- 
known problem, which is particularly severe in deep fields where 
both the background and the source number density are high. 
Therefore, using vtpdetect requires significant additional 
work to obtain a well-characterized list of extended sources. 

In this work we present a new end-to-end software 
by implementing the same VT+FOF method as used in 
Ebeling & Wiedenmann (1993) and solving these drawbacks. 
Most importantly, we add a self-contained deblending procedure 
and provide a well-characterized extended source list as the out- 
put. We preserve all advantages of vtpdetect, without intro- 
ducing any image binning or smoothing, or assuming any spe- 
cific source shape or size. The details of our algorithm are given 
in the following sections. 



2.2. Voronoi diagram construction 

The pattern of the Voronoi cells is also called a Voronoi diagram, 
which is composed of the edges of the cells (Voronoi edges). A 
simple algorithm for constructing a Voronoi diagram at the high- 
est computational speed is the sweep-line algorithm (Fortune 



1986) 3 . We implemented this algorithm and further developed 
it for better and faster performance. As a first step, we confined 
the Voronoi diagram, originally defined in the infinite plane, to 
the box domain of the image. Second, in addition to the clas- 
sical Voronoi diagram, we built a discrete Voronoi diagram by 
including all pixels enclosed by a Voronoi cell into a discrete 
Voronoi cell. If one pixel is equally distant from two or more 
Voronoi sites (i.e., the pixel center falls on an Voronoi edge), it 
is assigned randomly among the two or more cells. Finally, as a 
third step, we connected the pair of Voronoi sites on both sides 
of each Voronoi edge, building the so-called Delaunay diagram, 
which will constitute the structure on which FOF and our de- 
blending procedure are run. 

Based on the Voronoi diagram, we created an area map that 
contains the area of Voronoi cells associated to each occupied 
pixel (Voronoi site). We calculated the accurate area of the cell 
polygons rather than counting the pixels in the discrete Voronoi 
cells. Finally, we computed the photon density map by dividing 
the original image by the area map. An example of a Voronoi 
tessellation and its associated photon density map is shown in 
Figure 1. To identify source candidates, the photon density map 
must be compared with a background map, which is generated 
as described in the following section. 



2.3. Background map generation 

A photon in an X-ray image can be associated to a source 
or the background. The background is defined as the sum of 
all contributions not associated to astronomical sources, or as- 
sociated to some astronomical component that cannot be re- 
solved (like the Galactic diffuse emission). We estimated the 
average background flux adopting the same method as used 
in Ebeling & Wiedenmann (1993, see their section D). We as- 
sumed that the background photons in the exposure-corrected 
image, which is obtained by dividing the photon density image 
by the exposure map, are distributed according to the Poisson 
statistics, and that the fraction of filled pixels with the lowest 
photon density values are dominated by background. Fitting a 
simple analytical model to the distribution of the low photon- 
density pixels, an average background intensity was estimated. 
The final background map of the image in units of photon/ pixel 
is obtained by multiplying the mean background flux by the ex- 
posure map. In this way the variations across the field of view 
due to vignetting, edges of the chips, or missing columns and 
pixels were taken into account. 

We adopted a constant background for each exposure- 
corrected XRT image. This choice is reasonable for XRT images 
as we tested with simulations. However, it may not apply to im- 
ages with strong background variations. Another caveat is that 
this method works properly when the number of occupied pixels 
is dominated by background photons. Thus it cannot be applied 
to fields dominated by source emission (for example, very shal- 
low fields with very bright sources, or with extended sources that 
cover the entire FOV). For such fields a manual treatment of the 
background is preferable. We plan to introduce a more sophisti- 
cated background treatment in a future version of the code. 



3 The computing time scales as 0{n\ogri), where n is the number 
of occupied pixels. Using the sweep-line algorithm, the time spent on 
Voronoi construction is negligible compared to the subsequent steps. 
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Fig. 1: Left panel: original unbinned X-ray image with the Voronoi tessellation (green lines) obtained with the classical sweep-line 
algorithm. The value of each filled pixel is simply the number of photons. Right panel: photon density map obtained by dividing the 
original image by the area map. The intensity of each filled pixel is now proportional to the photon density. 



2.4. FOF detection and deblending 

In the FOF algorithm two occupied pixels whose distance is 
smaller than a chosen linking length are defined as friends. 
Starting from any occupied pixel, all its friends are associ- 
ated with it, and then friends of its friends, until no more new 
friends can be found. All connected pixels (or the corresponding 
Voronoi cells) constitute a source candidate. The FOF algorithm 
is iteratively applied to all involved pixels, until all of them are 
assigned to some source candidates. We ran FOF on the occupied 
pixels whose photon density was above the background level. 
This choice is very generous in the sense that several spurious 
detections may survive above the background level. However, 
we preferred to start with a large list of candidates and to refine 
it later, a procedure often adopted (see Broos et al. 2010). 

The direct application of the FOF algorithm results in 
a large number of blended sources, as mentioned also in 
Ebeling & Wiedenmann (1993). However, given the relatively 
low number density of sources in the X-ray sky 4 , most of the 
blended sources overlap only in their fainter outer regions , while 
their bright cores, corresponding to the peak of the PSF, remain 
isolated. In principle, the blended sources can be deblended by 
raising the threshold well above the background level. 



4 This clearly depends on the depth of the image. At the flux levels 
currently explored in the deepest survey (~ 10~ 17 erg cm -2 s _1 in the soft 
band in the CDFS, see Xue et al. 201 1), this value is about 2 x 10 4 deg" 2 
(Lehmer et al. 2012), and only thanks to the high angular resolution of 
Chandra this value is still far from the confusion limit. Nevertheless, 
most of the archival X-ray data available at present are, at best, one 
or two orders of magnitudes less sensitive than the CDFS, which will 
remain for a long time the deepest X-ray image of the sky. 



The threshold necessary to separate two blended sources 
varies according to the background, the distance between the 
sources, and the size of the sources themselves, which, for unre- 
solved sources, is determined only by the PSF profile. This step 
has been performed by Horner et al. (2008) on the vtpdet ect 
direct output. They ran vtpdet ect five times on each X-ray 
image using different thresholds of surface brightness, and then 
selected the best threshold for each field by visually inspect- 
ing the source photon distribution for each threshold value. This 
procedure is reliable but time-consuming, especially for a large 
survey that spans a wide range of exposure times and therefore 
surface brightness thresholds. As previously anticipated, a ma- 
jor improvement of our new software consists in including the 
auto-deblending procedure in the algorithm. 

As a first step, to flag blended sources, we created a number 
of source maps 5 by running FOF with thresholds correspond- 
ing to different photon densities, starting from the lowest (corre- 
sponding to the background level, map with index j = 1) to the 
highest value present in the photon density image, distributed 
on a logarithmic scale. Clearly, the source regions defined in the 
map j=l will be larger than the regions in the higher order 
maps. An example of this process is shown in Figure 2, where 
two blended sources are identified by FOF as one in panels 1 
to 4, and finally as two separate sources in panel 5. A source is 
tagged as blended whenever its source region includes more than 
one source in one of the higher order maps. 

At this point we had a list of sources virtually free from 
blending. Their positions were assigned by searching for the 
pixel whose 3x3 square island has the largest number of pho- 



5 The number of the maps can vary; a number around 10 guarantees 
robust results. 
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tons. Then we separated the blended sources among the most 
reasonable non-overlapping source regions, without assuming a 
source brightness model. This is obtained as follows. For each 
pair of blended sources, we can find a threshold above which 
they are detected without overlapping pixels. Sorting all occu- 
pied pixels by their photon densities, we identify the pixel with 
the lowest value whose inclusion causes the bridging of the two 
sources. This pixel is flagged as a bridge pixel and is removed 
from the original image. After a bridge pixel is found and re- 
moved, we again run FOF on all occupied pixels above local 
background. Unless the sources are detected separately by FOF, 
we repeat the above procedure to identify a new bridge pixel. 
The final set of bridge pixels is the minimum set of the faintest 
filled pixels we need to remove to separate the two sources. This 
process can be seen in Figure 2 going from panel 5 back to panel 
1 . As soon as the set of bridge pixels are found and removed, the 
regions of the two sources are built directly with FOF. At the last 
step, each one of the bridge pixels are re-assigned to one of the 
two sources, with a simple criterion based on the neighboring 
filled pixels: if the majority of them belong to a given source, 
the bridge pixel is assigned to this source. In the very few cases 
where this condition does not apply, the bridge pixel is assigned 
randomly. The final result is shown in panel 6 of Figure 2. 
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Fig. 2: Example of the deblending procedure applied to two un- 
resolved sources spaced by ~ 6 pixels. In panels 1-5 we show 
the pixels whose photon density is above a threshold that in- 
creases at each step. Bridge pixels (see text for the definition) 
are marked with green boxes and are identified backwards start- 
ing from panel 4 to 1 . Panel 6 shows the two sources identified 
separately, with the bridge pixels re-assigned according to a lo- 
cal criterion. 



There are some cases that require special treatments. The 
first one concerns blending of extended sources. In principle, 
blended extended sources cannot be separated without making 
assumptions on their surface brightness distributions. However, 
for most of the extended sources we are interested in, such as 
regular clusters with only one bright cool core, it is reasonable to 



assume that the surface brightness decreases monotonically from 
the core to the edge. In such cases, blended extended sources can 
be separated with the same method as described above. This step 
is left as an option. The second case is the bridging between ex- 
tended sources and background overdensity regions, which re- 
sults in branch-like shapes at the edge of extended sources. We 
managed to clean such cases by identifying the cut-vertex. A cut- 
vertex is similar to a bridge pixel, but it is defined as the single 
pixel whose removal allows one to separate the blended sources 
into two or more components in the image. Clearly, the pres- 
ence of a cut- vertex indicates a weak connection between two 
regions. Cutting off these weak connections allowed us to effi- 
ciently reduce the bridging of real extended sources with back- 
ground fluctuations. 

2.5. Reliability filter and source classification 

We applied a reliability criterion by comparing the total pho- 
tons C tota i inside the source region to the number of background 
photons Bkg expected in the same area, as obtained from the 
background map. We computed the S/N as C net / ^Bkg, where 
Cnet = C total ~ Bkg. This is not the traditional definition of the 
S/N, which is net counts over the square root of the total counts 
(source plus background) for Poisson statistics. However, this 
definition 6 is a good estimator of the probability of a source to be 
inconsistent with a Poissonian background fluctuation. In Figure 
3 we show the curves in the C net - Bkg plane that correspond to 
S/N = 3, 4 and to a Poissonian probability of 10" 3 and 10" 4 of to 
be a random fluctuations. A filter for a given reliability threshold 
can be applied by setting the S/N threshold to the required con- 
fidence level. In the remainder of the paper we set this threshold 
to S/N > 4, corresponding to a Poissonian probability of 10" 4 
of the signal inside the source region being random fluctuation. 
We verified a posteriori that with this condition our algorithm is 
able to detect extended sources down to a flux level much lower 
than that required to characterize them as extended. Therefore, 
this step does not significantly affect the final flux limit for ex- 
tended source detection. The filter was applied not only to the 
final source list, but also on the run to the temporary source list 
to avoid applying the extension criterion to unreliable sources 
and thus to reduce the computation time. 

To classify all sources according to their extent, we drew a 
circle centered on the position of each source after masking all 
surrounding detected sources. The radius was chosen to enclose 
60% ~ 70% energy (slightly larger than HEW/2 of the PSF). 
For the sources that survived the reliability filter, we produced a 
normalized cumulative profile in the core circle by sorting each 
pixel according to the distance from the source position, and 
summing all photons at shorter distance 7 . A reference profile for 
the instrument PSF was calculated in the same way. The two pro- 
files were then compared with a Kolmogorov-Smirnov (KS) test. 
Under the null hypothesis that they are drawn from the same dis- 
tribution, this procedure gives us the probability that the source 
image is not consistent with being unresolved in the core region. 
An example is shown in Figure 4. Based on this procedure, we 
can classify the source as unresolved or extended by setting a 
threshold on the null-hypothesis probability of the KS test. This 
threshold must be chosen a posteriori after extensive simulations 
to keep both the contamination and the completeness of the sam- 
ple under control. This is shown for Swiff-XRT in §3. 

6 See also http://heasarc.gsfc.nasa.gov/W3Browse/all/cxoxassist.html 

7 For an asymmetric PSF the distance of each pixel to the source 
position can be changed to the according PSF value of each pixel. 
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Pixels 



Fig. 4: Left panel: a simulated XRT image with five labeled sources, three AGNs (sources 0, 1, and 2) and two clusters (sources 
3 and 4). Right panel: the cumulative net count profile for the sources shown in the upper panel. The KS-test results are given for 
each source in the inset frame. The two clusters are clearly identified as extended (with a profile above the PSF and a high KS-test 
probability) and the three AGNs as unresolved sources (with a profile consistent with that of the PSF and a low KS-test probability). 
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Fig. 3: Curves in the net photons-background photons plane cor- 
responding to a probability p of 10" 3 or 10" 4 of being a ran- 
dom fluctuation (red lines), compared with those corresponding 
to S/N of 3 or 4 (blue lines). Curves for S/N = 4, p = 10" 4 
and S/N = 3, p = 10" 3 are shown with solid and dashed lines, 
respectively. 



The KS null-hypothesis probability also depends on the net 
detected photons and the background level. As we show in 
Section 3, one simple but effective way is to set the thresh- 
old that divides the unresolved from the extended source in 
the probability- S/N plane. The optimal threshold is determined 
through extensive simulations by computing the completeness 



(fraction of recovered extended sources) and the contamination 
(number of sources spuriously classified as extended) and choos- 
ing the best compromise. 

One may argue that a small region (the core circle) is not the 
best choice for classification, because extended sources are max- 
imally different from unresolved sources at large radii. However, 
because of the low S/N in the outer regions, it turns out that 
it is more efficient to focus on the core region where the S/N 
is maximized. However, this choice may cause the loss of ex- 
tended sources whenever a bright unresolved source in the cen- 
ter is present. 8 . This may occur because of a strong (typically 
unrelated) unresolved source is embedded in the extended emis- 
sion, or because of the presence of a bright cool core whose size 
is below the angular resolution. To overcome this problem, we 
applied a supplementary step to recover most of these sources. 
We considered all sources that are classified as unresolved but 
with KS null-hypothesis probability > 50%. For each of these 
sources, we normalized the local PSF to the total photon in- 
side the core region and calculated the radius where it reaches 
1/2 the local background level. Since the edges resulting from 
the FOF algorithm are set by the overall background level, the 
circles defined by this criterion should not be entirely encom- 
passed by the FOF source regions. If this happens, though, it 
may imply the presence of extended emission embedded in a 
larger background region (for whatever reason) and the source is 
tentatively classified with a flag as extended. We remark, how- 
ever, that this procedure may introduce several false candidates 
if the PSF shape is not completely under control. Eventually, 
these cases may be visually inspected for final inclusion in the 



8 This effect seems to be negligible in real X-ray images, as shown, 
for example, by the X-ray follow-up of optically and IR selected clus- 
ters, which indicate that the occurrence of strong X-ray unresolved 
sources associated with the extended emission from groups and clusters 
of galaxies is very rare (see Hicks et al. 2008; Bignamini et al. 2008) 
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extended source list. Finally, as a very last step, extended sources 
found within a larger, surrounding extended source are merged 
into it and considered part of the largest one. 

2.6. Source regions and aperture photometry 

The definition of a FOF detected source region depends on the 
linking length, which is the only parameter of FOF. The linking 
length is defined as the median distance between occupied pix- 
els in a pure background image, and therefore it depends on the 
background intensity only. An estimate of the linking length in 
units of pixel size is given by -\jBkg~ 1 , where Bkg is in units of 
photon/ pixel. Clearly the linking length depends on the position 
on the image, therefore we defined a linking length map accord- 
ing to the background map. Using this linking length, the bor- 
ders of each source are automatically defined as the loci where 
the photon density is equal to the background density. We recall 
that the photon density still includes both the source signal and 
the background, since no background was subtracted from the 
image. With this definition, statistically the entire signal associ- 
ated to the source is included in the FOF region. In other words, 
this choice of the linking length corresponds to a precise defini- 
tion of the source boundaries, and allows us to perform aperture 
photometry directly inside the FOF region. 

For the sources classified as extended, the regions obtained 
by FOF were used as photometry aperture. However, unresolved 
sources that are well-described by the PSF shape were treated 
differently. For each source classified as unresolved, we normal- 
ized the PSF to the total photon inside the core region and calcu- 
lated the radius where it reaches 1/2 the local background level. 
This defines the extraction radius for all unresolved sources. 
When two unresolved sources are closer than this radius, the 
overlapping pixels were assigned to the one with the brighter 
expected value in that position. A more sophisticated treatment 
including a source- fitting procedure is planned to be introduced 
in a future version of the algorithm. 

We calculated the total number of photons within each 
source region from the original image, and the total background 
photons within the same region from the background map. 
Aperture photometry for each source was simply computed as 
the total number of photons minus the number of background 
photons. After all these steps, we obtained a list of extended 
source candidates without any a priori assumptions on their in- 
trinsic size and properties. For each source, we provide position 
and area, together with its aperture photometry, and the normal- 
ized exposure time corresponding to the emission- weighted ex- 
posure in the source region. 



3. Application to Swift-XRT archival data 

One of our goals is to provide a flexible algorithm that can be 
used for extensive simulations adapted to specific cases. This is 
a relevant aspect since the performance of any algorithm changes 
significantly as a function of instrument characteristics, such 
as resolution and background level. In this section we describe 
set-up of simulations to obtain a characterization of a survey 
in terms of completeness and contamination. We describe the 
specific case of the X-ray Telescope (XRT) onboard the Swift 
satellite, which has been used mainly to follow up gamma-ray 
bursts (GRB), and which so far provided a serendipitous survey 
of about 70 deg 2 of the X-ray sky (see Tundo et al. 2012). The 
result of this section will be used in a forthcoming paper to char- 



acterize the final catalog (from about 400 deg 2 ) of SXCS (Liu et 
al. in preparation). 

3. 1. XRT image simulation 

Tundo et al. (2012) used a first chunk of the XRT archive as of 
April 2010 including 336 GRB follow-up fields with Galactic 
latitude \b\ > 20 deg, resulting in a catalog of 72 extended X- 
ray sources (catalog I). Several of these sources have already 
been identified as groups or clusters thanks to a cross-correlation 
with available optically data (mostly from the Sloan Digital Sky 
Survey) or have been confirmed as clusters thanks to a detailed 
X-ray spectral analysis (Moretti et al. in preparation). 

The selection method used in Tundo et al. (2012) was based 
on the wavdetect algorithm and a simple growth-curve char- 
acterization. Completeness and contamination levels were kept 
under control thanks to the relatively high detection threshold 
adopted in catalog I (corresponding to at least 100 net counts 
in the soft band). We intend to increase both the sensitivity and 
the sky coverage of SXCS by applying the algorithm described 
here to the full Swiff-XRT archive (updated to November 201 1), 
after removing all fields with Galactic latitude \b\ > 20 deg 
or originally targeted at groups or clusters. Thanks to a much 
larger number of fields (> 3000) and a lower flux limit than that 
adopted in Tundo et al. (2012), we expect a significantly larger 
number of extended sources (> 300). To push the detection limit 
of extended sources down to the faintest possible fluxes, we need 
to robustly assess the completeness and the contamination levels 
of the survey. 

To the aim of building simulations as close as possible to the 
actual survey, we generated one image for each field in the soft 
band using the same exposure time, exposure map, Galactic col- 
umn density, energy conversion factor (ECF), and background 
flux as in the real image (see Figure 5 for a comparison of a 
real image and a mock image of the same field). Point sources 
were randomly extracted from a distribution consistent with the 
logN-logS measured in deep Chandra fields (Rosati et al. 2002; 
Lehmer et al. 2012) and simulated down to fluxes about one or- 
der of magnitudes lower than the expected detection limit of 
XRT images. For each input source we converted the intrin- 
sic soft-band flux into photon rate, using the appropriate ECFs, 
which are computed as in Tundo et al. (2012), taking into ac- 
count the Galactic absorption. While computing the ECF, we 
assumed a spectral slope of T = 1.4, which is justified by 
the observed stacked spectrum of all unresolved sources in the 
flux range of interest, as shown in deep-field studies (see, e.g., 
Tozzi et al. 2001). Given the exposure time, we computed the 
expected photon rate from each simulated source. Then the im- 
age of each unresolved source was created by distributing these 
photons according to the PSF of XRT (Moretti et al. 2007). After 
randomly distributing the sources across the field, the exposure 
map was used to apply vignetting correction at a given position. 

It is admittedly a limitation of our simulations that we as- 
sumed no spatial correlation between sources. However, we ar- 
gue that adding a correlation among X-ray sources would not 
significantly affect the final results. In particular, our algorithm 
includes a treatment of unresolved sources embedded within 
extended ones (see the last paragraph in §2.5), which reduces 
the impact of the correlation between unresolved and extended 
sources. 

The background level of each field was measured from the 
actual SXCS image with the method described in Section 2.3. 
This background was added with a Poissonian distribution to the 
mock image. The added background already includes the con- 
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Fig. 5: Soft band (0.5-2.0 keV) image of the field 1 10059+514249 from the XRT archive (left), and a mock image simulated for the 
same field. 



tributions from sources below the detection limit of each image. 
Since the mock unresolved sources are simulated down to fluxes 
well below the detection limit in XRT images, the contribution 
of sub-threshold sources would be counted twice. To account for 
this effect, we consistently revised the value of the background 
downward. The background values were checked a posteriori by 
comparing the background of the final simulated images with 
that in the corresponding real images. 

The flux distribution of the input extended sources was taken 
from the logN-logS of groups and clusters measured in the 
ROSAT deep cluster survey (Rosati et al. 1998). We did not con- 
sider other extended sources such as radio jets or nearby resolved 
normal galaxies. As previously mentioned, the morphology of 
extended sources may affect detection at a significant level to- 
ward the faint end. To take into account the different morpholo- 
gies of groups and clusters, we used ten real images of relatively 
bright groups and clusters of galaxies obtained with the Chandra 
satellite, covering a wide range in ICM temperature (from 2 to 8 
keV). The source image was obtained by rescaling the original 
Chandra image to the number of net photons expected for XRT, 
and to the source size typically corresponding to the source flux, 
adopting the phenomenologically observed relation between X- 
ray size and flux. This is just an approximation adopted to avoid 
unrealistic sources with large size and low flux or vice versa, but 
the details or the intrinsic scatter of this relation do not affect 
the results of the simulations. This technique has already been 
used to investigate the evolution of cool core clusters at high-z 
in Santos et al. (2010). The expected net photons were computed 
as for the unresolved sources, adopting the ECF appropriate for 
a thermal spectrum, which, as shown in Tundo et al. (2012), has 



a weak dependence on the actual temperature and redshift of the 
source. Finally the mock cluster image was convolved with the 
XRT PSF. 

3.2. Detection efficiency and contamination 

According to the definition in §2.5, a core radius of 5 pixel sizes 
(corresponding to ~ 12 arcsec, enclosing ~ 60% of the total 
energy) was chosen. In Figure 6 we show the threshold in the 
KS-test probability-S/N plane that we used to classify a source 
as extended. In the same plane we show the difference between 
the number of extended and unresolved sources as a color-coded 
grid. This immediately shows that the threshold has been defined 
simply by the condition of having an equal number of spurious 
and extended source in a given position on the grid. In the fol- 
lowing we compute the fraction of recovered extended sources 
and the contamination level of our survey based on this simple 
choice. We stress that the final properties of the survey critically 
depends on this selection curve. One may want to change the se- 
lection curve to find the best compromise between completeness 
and contamination required for a given scientific goal: this can 
be done by applying different selection curves to the simulations. 

The completeness is defined as the fraction of recovered ex- 
tended sources in a given net photon bin. We have a 100% com- 
pleteness above ~ 200 input net soft photons, and a gentle de- 
crease down to 80% at ~ 60 photons. As shown in Figure 7, 
the completeness also depends on the source extent, with larger 
sources having higher identification probability. This effect of 
course depends on the cluster population we used, which may 
not reflect the actual cluster population at high-z (or faint fluxes) 
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Fig. 6: Selection curves in the KS-test probability-S/N plane 
(magenta solid line). The color-coded grid indicates the absolute 
value of the number of unresolved minus the number of extended 
sources. Regions dominated by unresolved sources are depicted 
in blue, while regions dominated by extended sources are shown 
in red. The selection curve has been defined as the loci of zeros 
on the grid. The cyan line refers to a KS-test probability of 50%. 



that is still largely unknown. In the following, when we consider 
the average completeness as a function of the source photometry, 
we account for the mix of surface brightness distributions used 
in the simulation. We acknowledge that a group/cluster popula- 
tion with different morphologies can give different completeness 
curves. 

The completeness curve can be used to set a net photon 
threshold to the final source list. For example, if one allows a 
minimum completeness level of 90%, all sources with more than 
80 net photons are included in the final catalog. The sharp limit 
in net photons immediately translates into a flux limit for each 
exposure time. For each field we computed a flux-limit map ob- 
tained as ECF(Nh) X &0/Expmap(t), where ECF depends on the 
Galactic column density, and Expmap(t) (the exposure map in 
units of effective time) includes the effect of vignetting. Then, 
the solid angle covered by the survey above a given flux was 
obtained by measuring the total solid angle where the flux-limit 
is lower than a given flux. In Figure 8 we show the resulting 
sky coverage £l(S ) of the total SXCS (Liu et al. in preparation) 
and compare it with the sky coverage of the 400 Square Degree 
ROSAT PSPC Galaxy Cluster Survey by Burenin et al. (2007) 
and with the sky coverage used in Tundo et al. (2012), which 
corresponds to a sharp limit of 100 net photons and a much 
smaller number of fields. 

The other key aspect to be kept under control is the con- 
tamination level. The cumulative number of unresolved sources 
spuriously classified as extended above a given photometry is 
shown in Figure 9. This value is rapidly increasing below 100 net 
counts. At around 80 net photons, the contamination is reaching 
a high level of about 40 sources in the entire survey. As shown 
in Figure 9, most of the contamination is due to deep fields (with 
exposure times > 10 5 s), which have a higher background and 
more blended sources. One option to reduce the contamination 
is to apply tighter constraints on the source classification. Of 
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Net input photons 

Fig. 7: Color-coded 2D completeness distribution in the 
photometry- size plane. For a given value of detected net pho- 
tons, larger sources are more easily identified (i.e., detected and 
characterized as extended). 




Flux erg cm 2 s 



Fig. 8: Sky coverage of the complete SXCS (Liu et al. in prepara- 
tion) compared with that of the first SXCS release (Tundo et al. 
2012) and with the 400sd survey (Burenin et al. 2007). 



course, this would affect the completeness, and the best com- 
promise must be fine-tuned with several trials. Another option is 
to proceed with a visual inspection of all extended source candi- 
dates. 

Visual inspection has already been used in Tundo et al. 
(2012), and it is often used in XMM-Newton cluster surveys as 
well. It consists of a careful image inspection of each extended 
source candidate to reject all sources that appear to be domi- 
nated by the contribution of unresolved sources. Since most of 
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the spurious candidates are found in high-exposure fields, this 
procedure is performed only on fields with exposure times > 10 5 
s. As shown in Figure 9 by comparing the solid and dashed lines, 
this procedure reduces the number of spurious sources with more 
than ~ 80 photons to an acceptable level: only ~ ten spurious 
sources are now expected in the entire survey (about 400 deg 2 ), 
as opposed to 40 before the visual inspection. In our simulations, 
only one extended source with more than 80 net photons has 
been erroneously discarded with the visual inspection. This re- 
sult shows that X-ray astronomy would greatly benefit by visual 
analysis, as already successfully done for optical extragalactic 
astronomy in the form of "crowd sourcing" experiments (see, 
e.g., Lintottetal. 2008). 
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Fig. 9: Cumulative number (top) and fraction (bottom) of sources 
spuriously classified as extended as a function of the measured 
net photons, before (blue) and after (red) visual inspection. 
Results of fields with exposure time < 10 5 s and > 10 5 s are plot- 
ted separately, showing that most of the contamination comes 
from deep fields. 



We also verified the accuracy of our extended source pho- 
tometry. In Figure 10 we show the input versus the recovered 
net photons in the soft band for all sources identified as ex- 
tended. Spurious extended sources (marked as blue points) typ- 
ically have overestimated fluxes, because they are often due to 
blending of unresolved sources. For all properly characterized 
extended sources the photometry agrees very well with the ac- 
tual fluxes and no significant bias is observed. The cumulative 
number counts for sources detected with more than 80 net pho- 
tons is finally computed as 



N(>S) = Z Si >sC7 l /n(S i ), 



(D 



sponding to Si (see Tundo et al. 2012). Finally, Q is the average 
completeness factor corresponding to a given source photome- 
try, which depends on the net detected photons of the I th source 
and it is estimated on the basis of the mix of surface brightness 
distribution assumed in the calibrating simulations (see Figure 
7). In this way, each source is weighted with a factor inversely 
proportional to the survey completeness. 

This treatment is correct only if we can ignore the effects 
of the statistical errors on the flux measurements. In reality, 
flux errors, due to the Poissonian noise on the aperture pho- 
tometry of each source, introduce a Malmquist bias. In other 
words, the measured cumulative number counts N(> S) are bi- 
ased high by an amount roughly estimated as AN(> S ) ~ axN(> 
S)ACts/Cts, where a is the slope of the logN-logS, and Cts is 
the typical number of measured photons at the survey limit. If 
we assume a low background, as in the case of XRT, we can as- 
sume ACts ~ yjCts. For Cts ~ 80 and N ~ 300 as in our case, 
we estimate AN ~ 30, which implies about 35 sources above 
the detection limit due to the Malmquist bias. This is a crude 
estimate, but the Malmquist bias is consistently included in our 
simulations, since it can be directly computed from the ratio of 
the recovered and input sources as a function of the measured to- 
tal net photons. Once the Malmquist bias was taken into account, 
we checked that we recovered the input model of the logN-logS 
with good accuracy. 




where S is the total soft flux, S [ is the soft flux within the extrac- 
tion radius of the I th source, and Q(S;) is the sky coverage corre- 



Input net count 

Fig. 10: Net recovered photons vs input photons for the extended 
sources detected in the simulation (red points). Spurious ex- 
tended sources are also shown (blue points). The best-fit for 
sources with more than 80 net photons is C out = 0.996 x Cm + 
0.383. 



4. Improving the angular resolution 

Angular resolution is a key factor for detecting faint sources. A 
high angular resolution enables the identification of unresolved 
sources even with an extremely low number of photons (about 5 
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for an angular resolution of 1" - 90% encircled energy width, so 
far achieved only by Chandra at the aimpoint, see Broos et al. 
2010). Angular resolution is clearly also crucial to distinguish 
between extended and unresolved X-ray sources, and to remove 
the contribution of faint unresolved sources, which otherwise 
would be confused with the extended emission itself. To study 
the impact of angular resolution on X-ray surveys and in par- 
ticular in the identification of extended sources, we again ran 
the entire set of simulations after improving the PSF by a fac- 
tor of two. This corresponds to an imaginary instrument with 
the same properties as Swift-XRT, including pixel size, field of 
view, and background level, but its PSF HEW is only 9 arcsec, 
half of the original HEW of XRT (Moretti et al. 2007). With this 
sharper PSF, the core radius was chosen to be 3 pixel sizes (~ 7 
arcsec, enclosing ~ 65% of the total source energy of an unre- 
solved source). A classification threshold curve was chosen with 
the same criterion. As shown in Figure 11, the completeness 
of the simulation is significantly improved, while the contami- 
nation is reduced. With this enhanced resolution, we can iden- 
tify extended sources down to a much lower flux limit at similar 
completeness and contamination levels. We did not include the 
visual inspection here, but instead compared the direct outcome 
of the source detection algorithm. As a rule of thumb, we find 
that by improving the angular resolution by a factor of two we 
can identify extended sources down to a flux limit about 60% 
lower for the same completeness and contamination levels. We 
also expect that visual inspection would improve the final results 
significantly. 

This exercise shows that angular resolution is critical for 
the performance of an X-ray survey mission. While a power- 
ful, small FOV, pointing telescope should maximize angular res- 
olution at the aimpoint as in the classical Wolter type mirrors, 
a survey telescope should maximize the area-weighted angular 
resolution on the entire FOV. The only mission designed accord- 
ing to this criterion is the proposed Wide Field X-ray Telescope 
(WFXT, Murray & WFXT Team 2010; Rosati et al. 201 1). The 
advantage with respect to a Wolter type-I mirror configuration 
is clearly shown by the comparison of the 1 deg 2 Chandra 
COSMOS image, worth of the 1.8Ms (Elvis et al. 2009), to 
the simulated WFXT COSMOS field, obtained with only 13ks 
WFXT exposure (Tozzi et al. 2011). The simulated WFXT im- 
age has an angular resolution lower by only a factor of 2 than the 
average resolution of the Chandra mosaic, in front of a ~ 10 time 
larger HEW of the WFXT PSF with respect to Chandra. This is 
because using a Wolter type-I telescope in survey mode results 
in an average resolution much poorer than the nominal aimpoint 
resolution, due to the rapid degradation of the angular resolution 
at large off-axis angle, where most of the FOV resides. 

5. Conclusions 

The combined use of VT and FOF can provide an efficient al- 
gorithm for detecting faint extended sources in X-ray images. 
Although such an algorithm is available to the scientific commu- 
nity for many years (Ebeling & Wiedenmann 1993), it has not 
been used widely for searching extended sources, possibly be- 
cause of the complex implementation of the method. In this work 
we present an updated implementation of VT and FOF, plus an 
automated deblending procedure, in a software (EXSdetect) for 
the identification of extended source detection in X-ray images. 
The aim is to provide a user- friendly, end-to-end algorithm that 
can be used to exploit present X-ray data archives and to explore 
the performance of future X-ray missions. Among the most rel- 
evant properties of EXSdetect are the following: 
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Fig. 1 1 : Distributions of differential recovered fraction (upper 
panel) and cumulative contamination fraction (lower panel) in 
our simulations, obtained with standard XRT PSF (HEW =18 
arcsec, blue line) and with a sharper PSF (HEW = 9 arcsec, red 
line). 



- we include source detection, classification and photome- 
try, which are usually performed independently, in a single 
stand-alone algorithm. We also include a self-consistent de- 
blending procedure that efficiently reduces the number of 
blended sources; 

- thanks to our python implementation of the sweep-line algo- 
rithm, the adopted Voronoi construction scheme is fast also 
on large images; 

- only the PSF of the instrument across the FOV is necessary 
to run the algorithm, and no a priori assumption on the shape 
of the extended sources is needed; 

- the reasonably short computation time allows one to run ex- 
tensive simulations to optimize the internal parameters such 
as the selection curve for detecting and characterizing ex- 
tended sources, which depend on the survey properties. 

We also tested our algorithm on extensive simulations run 
for the SXCS survey. While the survey presented in Tundo et al. 
(2012) consists of about 40 deg 2 , here we considered a substan- 
tial extension of the survey, with a huge increase of the sky cov- 
erage particularly at bright fluxes. We find that, with our algo- 
rithm, we can detect and characterize extended sources with total 
net photons as low as 80 (in the soft band) with a completeness 
higher than 90%. The contamination amounts to a few tens of 
sources, but can be drastically reduced by visual inspection. In 
the simulations we recovered the input logN-logS with great ac- 
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curacy down to a flux of 10" 14 erg cm -2 s" 1 . This algorithm is 
currently being applied to the real SXCS data and the final cat- 
alog will be presented in Liu et al. (in preparation). The simple 
exercise of improving the angular resolution by a factor of two 
shows that the corresponding sensitivity to extended source de- 
tection improves by a factor -1.6. 

To summarize, EXSdetect is a new tool to exploit the 
huge X-ray archives from existing X-ray facilities like Chandra, 
XMM-Newton, Swift, and Suzaku. At the same time it is also 
very useful to explore the capability of future X-ray facilities. 
The Python code EXSdetect is available on the SXCS web- 
site (http : / /adlibitum .oats.inaf.it/sxcs) and it 
is open to continuous refinements and updates. 
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